function chi2 = minimization_function_v4(x, transitions_matrix, force_rate, force_distribution, detect_fluxes, variables_indexes,  initial_rates, chi2_minimization)

disp(num2str(x))
% optimization function returns the absolute difference between theoretical
% and experimental distributions
% x is the vector with the rates parameters: unstrained rates and force suceptibility   
    chi2 = 0;
    
        
    for m = 1:length(force_rate)
        initial_conditions = zeros(size(transitions_matrix, 1), 1);
        initial_conditions(1) = 1;
        non_zero_index = find(force_distribution(:,m+1)>0);
        time_interval =[0 (force_distribution(non_zero_index(end),1)+3)/force_rate(m)];
%         time_interval =[0 force_distribution(end,1)/force_rate(m)];
        %     time_interval =0 : 0.002 : (2 * max(force_distribution(:, 1)) / force_rate);

        options = odeset('RelTol', 1e-5, 'AbsTol', 1e-6, 'MaxStep', 0.1*5/force_rate(m), 'Refine', 10);
        %df = @(t, y)master_function_v4(t, y, x, transitions_matrix, force_rate(m), variables_indexes,  initial_rates);
        df = @(t, y)master_function_v4_mex(t, y, x, transitions_matrix, force_rate(m), variables_indexes,  initial_rates);


    %     disp('Integrating ...');
    
        %[T, Y] = ode45(df, time_interval, initial_conditions, options);
        [T, Y] = ode15s(df, time_interval, initial_conditions, options);

        N = length(initial_rates)/2;
        flux = zeros(length(T), N);
        for k = 1 : (length(T) - 1)
            s = 1;
            for i = 1 : length(transitions_matrix)
                for j = 1 : length(transitions_matrix)
                    if(transitions_matrix(i, j) > 0)
                        
                        index = find(variables_indexes==s); 
                        if(sum(index)>0)
                            k_0 = x(index);
                        else
                            k_0 = initial_rates(s);
                        end

                        index = find(variables_indexes==(s+N)); 
                        if(sum(index)>0)
                            x_0 = x(index);
                        else
                            x_0 = initial_rates(s+N);
                        end
                
                        rate = k_0 * exp(T(k) * force_rate(m) * x_0 / 4.1);
                       
                        if(rate >= 5000)
                            rate = 5000;
                        end
                        flux(k, s) = rate * Y(k, j);
                        s = s + 1;
                    end
                end
            end
        end

        if(length(detect_fluxes) > 1)
            time_distribution = sum(flux(:, detect_fluxes)');
        else
            time_distribution = flux(:, detect_fluxes)';
        end

%         time_distribution(time_distribution<0.002) = 0;
        time_distribution = time_distribution';
        normalizing_constant = cumtrapz(T(1:end), time_distribution(1:end));
    %     normalizing_constant(end)
        time_distribution = time_distribution/normalizing_constant(end);

        bin_size = force_distribution(2, 1) -  force_distribution(1, 1);
        histogram_start =  force_distribution(1, 1) - bin_size / 2;
        histogram_end =  force_distribution(length(force_distribution), 1) + bin_size / 2;

        events_number = sum(force_distribution(:, m+1));
        force_distribution(:, m+1) = force_distribution(:, m+1) / events_number;

        number_of_columns = (histogram_end - histogram_start) / bin_size;
        norm_coeff = events_number / (number_of_columns - length(x) - 1);
        theoretical_histogram = zeros((histogram_end - histogram_start) / bin_size, 1);
       difference = zeros((histogram_end - histogram_start) / bin_size, 1);
        for i = 1 : number_of_columns
            temp = [];
            t_start = (histogram_start + bin_size * (i - 1)) / force_rate(m);
            indexes = intersect(find(T >= t_start), find(T <= t_start + bin_size / force_rate(m)));
            temp = cumtrapz(T(indexes), time_distribution(indexes));
            theoretical_histogram(i) = temp(end);
            if(chi2_minimization)
                if( force_distribution(i, m+1) > 5/events_number)
                    difference(i) = norm_coeff * (theoretical_histogram(i) - force_distribution(i, m+1))^2/theoretical_histogram(i);
                    chi2 = chi2 + abs(difference(i));
                end
            else
                difference(i) = abs(theoretical_histogram(i) - force_distribution(i, m+1));
                chi2 = chi2 +difference(i);
            end
        end
 
%         cumtrapz(force_rate(m)*T, time_distribution/force_rate(m))
% chi2_1 = force_rate(m)*T;
% chi2_2 = bin_size*time_distribution/force_rate(m);
% chi2 = [chi2_1 chi2_2];
%         theoretical_histogram
%          force_distribution(:, m+1)
    %     disp(['delta =' num2str(y)]);

        figure(2*m-1)
        cla
%         plot(force_distribution(:,1), difference);
%         hold on
       plot(force_rate(m)*T, flux);

        figure(2*m)
        cla
        plot(force_distribution(:,1), theoretical_histogram, 'b');
        hold on
        plot(force_distribution(:,1), force_distribution(:,m+1), 'k');
        hold off

    end       
    
%     save_data = [force_distribution(:,1) force_distribution(:,m+1) theoretical_histogram];
%     save('c:\temp\1.txt', 'save_data', '-ascii', '-double', '-tabs')
    
%     figure(1)
%     hold on
%     ColorOrder = [0 0 1; 0 0.75 0; 1 0 0; 1 0.75 0; 1 0 0.75; 0.75 1 0; 0 1 0.75; 0.75 0 1; 0 0.75 1; 0 0.75 0.75; 0.75 0 0.75; 0.75 0.75 0];
%     for i = 1:size(Y,2)
%         color = ColorOrder(rem(i,12)+1,:);
%         plot(T, Y(:,i), '-', 'Color',color);
%     end
%     
%     figure(2)
%     hold on
%     ColorOrder = [0 0 1; 0 0.75 0; 1 0 0; 1 0.75 0; 1 0 0.75; 0.75 1 0; 0 1 0.75; 0.75 0 1; 0 0.75 1; 0 0.75 0.75; 0.75 0 0.75; 0.75 0.75 0];
%     for i = 1:size(flux,2)
%         color = ColorOrder(rem(i,12)+1,:);
%         plot(T, cumtrapz(T,flux(:,i)), '-', 'Color',color);
%     end
end

